# FDR, P-value, Q-value calculation
# 2014/01/30

# Usage
# Rscript Statistics.R [significant region file]

File <- commandArgs(trailingOnly =T)[1]
data <- read.table(File, header=F, sep="\t")

# file format
#chr start end MAXpos treat control SUM_treat SUM_control


### p-value
p <- pbinom(round(data$V7), round(data$V7 + data$V8), 0.5, lower.tail=F)



### Q-value
library(qvalue)
qobj <- qvalue(p)
q <- qobj$qvalues


### fold-change
fc <- data$V7 / data$V8


### output
cat("#chr", "start", "end", "MAX_pos", "treatment", "control", "readPerBp_treatment", "readPerBp_control", 
	"fold-change", "P-value", "Q-value\n", sep="\t", file=File, append=FALSE)
output <- data.frame(data, foldChange=fc, Pvalue=p, Qvalue=q)
write.table(output, file=File, append=TRUE, quote=F, sep="\t", row.names=F, col.names=F)

